QIIME 2用户文档. 7差异丰度分析gneiss(2018.11)

前情提要

QIIME 2用户文档. 7差异丰度分析gneiss

Differential abundance analysis with gneiss

原文地址： https://docs.qiime2.org/2018.11/tutorials/gneiss/

此实例需要一些基础知识，要求完成本系列文章前两篇内容：《1简介和安装》和《4人体各部位微生物组分析》。

在本教程中，您将学习如何使用gneiss中的balances来进行差异丰度分析。我们将关注的主要问题是如何以组成连贯的方式识别差异丰富的分类群。

详者注：balances批比值中的分子，对应中文中的平衡、均衡、天平或余额都有相似，但又不完全一致。目前没有标准翻译，文中还是使用原词 balances。

组成性指的是处理比例的问题。为了考虑测序深度的差异，微生物的丰度通常被标准化为比例（例如，相对丰度）。正因为如此，准确地推断出哪些微生物正在发生变化就变得具有挑战性——因为比例加和为1，单个微生物的变化也会改变其余微生物的比例。

考虑以下示例：

图1. 相对丰度变化示意图

在左边，我们有十个物种的真实丰度，第一个物种在时间点1和时间点2之间加倍。当我们将这些按比例归一化时，似乎所有的物种在这两个时间点之间都发生了变化。单从比例上看，我们永远不会意识到这个问题，实际上我们不能仅仅根据比例精确地确定哪些物种正在发生变化。

虽然我们不能确切地解决鉴别不同物种的问题，但我们可以放宽这个问题，并询问微生物的哪些分区/部分(partitions)正在改变。在上面的例子中，如果我们计算第一种和第二种之间的比率，对于原始丰度和比例，这个比率在时间点1是1:1，在时间点2是2:1。这是balances试图解决的问题。我们可以把重点放在分类群（或分类组）之间的比率上，而不是放在单个分类群上，因为这些比率是由所观察物种的真实丰度和所观察的比例构成的。我们通常对这些比率进行对数转换，以优化可视化效果（“对数比”）。计算多个物种的balances（或ratios）的概念可以扩展到树，如下面的示例所示。

图2. 树形图展示原始值、相对丰度、比例间的关系

在左边，我们定义一棵树，其中每个树枝尖对应一个分类单元，下面是第一个样本中每个分类单元的比例。内部节点（即balances）定义了底下分类群之间的对数比率。右边是同一棵树，下面是不同样本中每个分类群的比例。只有一个分类群数量变化，正如我们之前所观察到的，所有分类群的比例都会改变，但是看看balances，只有包含紫色分类群的平衡才会改变。在这种情况下，balances不会改变，因为它只考虑红色与分类群之间的比率。通过观察balances而不是比例，我们可以通过限制观察只关注给定balances内的分类群来消除一些差异。

这里突出的问题是，我们如何构造balances树来控制这种变异，并识别出有趣的分类群差异丰富的分区？在gneiss中，这有三种主要的方法：

• 相关聚类(Correlation clustering)。如果我们没有关于如何将生物体聚集的相关先验信息，我们可以根据它们彼此共生的频率来将生物体分组在一起。这可以在correlation-clustering命令中获得，并为ilr-hierarchical创建树输入文件。

• 梯度聚类(Gradient clustering)。使用元数据类别对在相似样本类型中发现的分类群进行聚类。例如，如果我们要评估pH是否是驱动因子，我们可以根据观察到的分类群的pH进行聚类，并观察低pH生物与高pH生物的比例是否随着pH的变化而变化。这可以在gradient-clustering命令中获得，并为ilr-hierarchical创建树输入文件。

• 系统发育分析(Phylogenetic analysis)。也可以使用gneiss以外创建的系统发育树（例如，rooted-tree.qza）。在这种情况下，您可以使用您的系统发育树作为ilr-phylogenetic的输入。

一旦我们有了一棵树，我们就可以使用下面的等式来计算balances：

bi=rsr+s_logg(xr)g(xs)

其中i表示树中的第i个内部节点，g(x)表示集合x内的几何平均值，x r表示balances分子中的分类群丰富度集合，x s表示balances分母中的分类群丰富度集合，r和s分别表示x r和x s内的分类群数量。

在计算出balances后，可以进行标准统计过程，如方差分析和线性回归。我们将使用慢性疲劳综合症(Chronic Fatigue Syndrome )数据集演示运行这些过程。

创建balances

Creating balances

在Giloteaux等人(2016)发表的慢性疲劳综合征数据集中，有87人，48名患者和39名健康对照组。本教程使用的数据是在Illumina MiSeq上使用地球微生物组项目高变区4 (V4) 16S rRNA基因测序方法产出的。

对于上文提到了两种常用安装方法，我们每次在分析数据前，需要打开工作环境，根据情况选择对应的打开方式。

启动工作环境并创建工作目录

# 创建qiime2学习目录并进入 mkdir -p qiime2 cd qiime2 # Miniconda安装的请运行如下命令加载工作环境 source activate qiime2-2018.11 # 如果是docker安装的请运行如下命令，默认加载当前目录至/data目录 # docker run --rm -v $(pwd):/data --name=qiime -it  qiime2/core:2018.11 # 创建本节学习目录 mkdir qiime2-chronic-fatigue-syndrome-tutorial cd qiime2-chronic-fatigue-syndrome-tutorial

实验数据下载

注意：QIIME 2 官方测试数据部分保存在Google服务器上，国内下载比较困难。可使用代理服务器(如蓝灯)下载，或公众号后台回复”qiime2”获取测试数据批量下载链接，你还可以跳过以后的wget步骤。

下载来源Google文档的实验设计

wget \  -O "sample-metadata.tsv" \  "https://data.qiime2.org/2018.11/tutorials/gneiss/sample-metadata.tsv"

下载特征表和物种注释

wget \  -O "table.qza" \  "https://data.qiime2.org/2018.11/tutorials/gneiss/table.qza" wget \  -O "taxa.qza" \  "https://data.qiime2.org/2018.11/tutorials/gneiss/taxa.qza"

首先，我们将定义我们想要构建balances的微生物的分区。同样，有多种可能的方法来构建一个树（即层级结构），它定义了我们要为其构建balances的微生物balances的分区。我们将在这个数据集中展示相关聚类( correlation-clustering)和梯度聚类(gradient-clustering)。

请注意，我们将运行的差异丰度技术将使用对数比转换(log ratio transforms)。由于不允许取零的对数，下面的两种聚类方法都包含一个默认的伪计数参数。这将用1替换表中的所有零，这样我们就可以在转换后的表上应用对数。
输入表是原始计数表（FeatureTable[Frequency]）。

选项1：相关性聚类

Option 1: Correlation-clustering

这个选项应该是您的默认选项。我们将通过Ward的分层聚类来采用无监督聚类，以获得Principal Balances。从本质上讲，这将使用Ward层次聚类来定义通常彼此共存的微生物的分区，这是由以下度量标准定义的。

D（x，y）＝V[Ln(x/y)]

其中x和y代表两种微生物在所有样品中的比例。如果两种微生物高度相关，那么这个数量将缩小到接近零。然后，Ward层次集群将使用这个距离度量来迭代地将相互关联的微生物群聚集在一起。最后，我们得到的树将突出显示高层结构，并识别数据中的任何块。

qiime gneiss correlation-clustering \  --i-table table.qza \  --o-clustering hierarchy.qza

输出对象:

• table.qza: 特征表

• taxa.qza: 物种注释

• hierarchy.qza: 层级结构

选项2：梯度聚类

Option 2: Gradient-clustering

相关聚类的另一种选择是基于数字元数据类别创建树。通过梯度聚类，我们可以对出现在元数据类别类似范围内的分类群进行分组。在本例中，我们将使用元数据类别年龄创建树（层次结构）。请注意，元数据类别不能缺少变量，并且必须是数字。

qiime gneiss gradient-clustering \  --i-table table.qza \  --m-gradient-file sample-metadata.tsv \  --m-gradient-column Age \  --o-clustering gradient-hierarchy.qza

输出对象:

• gradient-hierarchy.qza: 梯度层级

下游分析的一个重要考虑因素是过度拟合问题。当使用梯度聚类时，您正在创建一个树来最好地突出所选元数据类别中的成分差异，并且可能会得到假阳性结果。小心使用梯度聚类。

用平衡建立线性模型

Building linear models using balances

现在我们有了一个定义分区的树，我们可以执行等轴测对数比率（ILR）转换。ILR转换计算树中每个节点上组之间的对数比率。

qiime gneiss ilr-hierarchical \  --i-table table.qza \  --i-tree hierarchy.qza \  --o-balances balances.qza

输出对象:

• balances.qza 对数比结果

既然我们有了树中每个节点的对数比率，我们就可以对balances进行线性回归。我们将要运行的线性回归称为多变量响应线性回归(multivariate response linear regression)，其归根结底是对每个balances分别执行线性回归。

我们可以用它来尝试将微生物丰度与环境变量联系起来。与标准单变量回归相比，运行此模型有多方面的优势，因为它可以避免与过度拟合相关的许多问题，并且可以让我们根据环境参数了解群落范围内的扰动。

由于微生物丰度可以直接映射到balances上，我们可以直接对balances进行多元回归。我们将要建立的模型表示如下

y→=β0→+βSubject→Xsubject→+βsex→Xsex→+βage→XAge→+βsCD14ugml→XsCD14ugml→+βLBPugml→XLBPugml→

其中y是代表待预测的balances矩阵，βi→表示协变量i的系数向量，Xi→表示协变量i的测量。

记住，方差分析(ANOVA)是线性回归的一种特殊情况——方差分析可以解决的每个问题都可以重新表述为线性回归。详情请参阅本帖。所以我们仍然可以用这种方法回答同样的差异丰度问题，我们可以控制多个不同的潜在混淆变量和交互效应以获得更精确的答案。

为了更好的理解以下公式，我们先看一下实验设计的样式

less -S sample-metadata.tsv#SampleID       Sample_Name_s   BarcodeSequence LinkerPrimerSequence    Subject Sex     Age     Pittsburgh      Bell    BMI     sCD14ugml       LBPugml LPSpgml IFABPpgml       Physical_functioning    Role_physical   Role_emotional  Energy_fatigue  Emotional_well_being    Social_functioning      Pain    General_health  Description ERR1331797      LR17    AAGGGACAAGTG    TATGGTAATTGT    Patient Female  67      2       40.0    32.89   1.65    15.01   279.3   145.0   40.0    0.0     0.0     35.0    52.0    25.0    23.0    30.0     ERR1331823      IC19    TCCACCCTCTAT    TATGGTAATTGT    Control Female  43                      26.15   1.36    18.76   100.54  184.1                                                                   ERR1331841      LR49    CTGTAGCTTGGC    TATGGTAATTGT    Control Female  60                      30.89   1.23    12.43   97.36   234.4

线性回归分析：

qiime gneiss ols-regression \  --p-formula "Subject+Sex+Age+BMI+sCD14ugml+LBPugml+LPSpgml" \  --i-table balances.qza \  --i-tree hierarchy.qza \  --m-metadata-file sample-metadata.tsv \  --o-visualization regression_summary.qzv

输出对象:

• regression_summary.qzv 回归摘要

图5. 回归结果总结。说明见下方段落。

现在我们总结了回归模型。具体来说，我们想看看哪些协变量(covariates)对模型影响最大，哪些balances有意义，以及有多少潜在的过度拟合正在发生。

在回归摘要中，有几点需要注意。摘要中有一个R2(Rsquared)，它给出了有关回归模型解释群落中有多少方差的信息。从我们看到的，回归可以解释大约11%的群落结构变异(Rsquared:    0.1108)。这是典型的人类肠道微生物群，因为由于遗传、饮食、环境等因素，有非常高的混杂变异。

为了评价一个单一协变量的解释模型，采用了一个遗漏变量的方法。遗漏一个变量，计算r2的变化(R2diff)。变化越大，协变量的作用越强。在这种情况下，课题(Subject)是最大的解释因素，解释了2.41%的变化。

为了确保我们不过度拟合，进行了10倍交叉验证。这将把数据分成10个分区，在其中9个分区上构建模型，并使用剩余的分区来度量预测的准确性。此过程重复10次，每轮交叉验证一次。每轮交叉验证报告模型内误差（within model error, mse）、R2和预测精度（R2 and the prediction accuracy, pred_mse）。在这里，预测精度低于模型内误差，表明不会发生过拟合。

接下来，我们采用热图可视化每个balances的所有协变量的显著性水平(p值)。热图的列代表balances，热图的行代表协变量。热图由p值的负对数着色，突出潜在的显著p值。使用鼠标悬停功能可以获得特定的系数值及其对应的p值，并启用缩放可以探索感兴趣的协变量和balances。

接下来是预测(prediction)和残差(residual)图。这里，只绘制了前两个balances，模型中的预测残差被投影到这两个balances上。从这些图中我们可以看出，预测点与原始群落位于同一区域内。然而，我们可以看到，残差与预测的方差大致相同。这有点让人不安——但请注意，我们只能解释10%的群落差异，所以这些计算很正常。

可视化树中的分支长度也可以通过模型中可解释的平方和(sum of squares)进行标准化。最长的分支长度对应于最有用的balances。这可以让我们对模型中最重要的balances进行高层次的概述。从这个图和上面的热图中，我们可以看到balances y0很重要。这些balances不仅具有很小的p值（p<0.05），能用于区分受试者，而且它们在树形图中具有最大的分支长度。这表明微生物的这种划分可以区分CFS患者和对照组。

我们可以在一个热图上可视化这些balances，以了解它们代表的分类群。默认情况下，特征表中的值是按对数标准化的，意味着样本值以零为中心。

qiime gneiss dendrogram-heatmap \  --i-table table.qza \  --i-tree hierarchy.qza \  --m-metadata-file sample-metadata.tsv \  --m-metadata-column Subject \  --p-color-map seismic \  --o-visualization heatmap.qzv

输出可视化结果:

• heatmap.qzv 热图

如图例所示，每个balances的分子(numerators)以浅红色突出显示，而分母以深红色突出显示。从这里我们可以看到，y0的分母与y0的分子相比，几乎没有OTU。分母中的这些分类群可能重要，所以让我们研究一下用balance_taxonomy构成balance的分类。

具体来说，我们将绘制一个箱线图，并确定能够解释对照组和患者组之间差异的分类群。

qiime gneiss balance-taxonomy \  --i-table table.qza \  --i-tree hierarchy.qza \  --i-taxonomy taxa.qza \  --p-taxa-level 2 \  --p-balance-name 'y0' \  --m-metadata-file sample-metadata.tsv \  --m-metadata-column Subject \  --o-visualization y0_taxa_summary.qzv

输出可视化结果:

• y0_taxa_summary.qzv 热图

在这种特殊情况下，与对照组相比，患者组的对数比率更低。实质上，这意味着与患者组相比，健康对照组y0分子中的分类群平均比y0分母中的分类群更丰富。

记住，根据本教程开头给出的示例，不可能推断出给定样本中微生物的绝对变化。Balances将无法提供这类答案，但它可以限制可能出现情况的数量。具体来说，可能发生了以下五种情况之一。

1. 患者组与健康对照组y0分子的平均分类数增加。

2. 患者组与健康对照组间y0分母的分类平均下降；

3. 上述两种情况的结合；

4. y0分子和y0分母的类群丰度均增加，但y0分子的类群丰度比y0分子增加更多。

5. y0分子和y0分母的类群丰度均下降，而y0分子比y0分母的类群丰度相对增加较多。

为了进一步缩小这些假设，需要生物先验知识或实验验证。

Reference

Bolyen E, Rideout JR, Dillon MR, Bokulich NA, Abnet C, Al-Ghalith GA, Alexander H, Alm EJ, Arumugam M, Asnicar F, Bai Y, Bisanz JE, Bittinger K, Brejnrod A, Brislawn CJ, Brown CT, Callahan BJ, Caraballo-Rodríguez AM, Chase J, Cope E, Da Silva R, Dorrestein PC, Douglas GM, Durall DM, Duvallet C, Edwardson CF, Ernst M, Estaki M, Fouquier J, Gauglitz JM, Gibson DL, Gonzalez A, Gorlick K, Guo J, Hillmann B, Holmes S, Holste H, Huttenhower C, Huttley G, Janssen S, Jarmusch AK, Jiang L, Kaehler B, Kang KB, Keefe CR, Keim P, Kelley ST, Knights D, Koester I, Kosciolek T, Kreps J, Langille MG, Lee J, Ley R, Liu Y, Loftfield E, Lozupone C, Maher M, Marotz C, Martin BD, McDonald D, McIver LJ, Melnik AV, Metcalf JL, Morgan SC, Morton J, Naimey AT, Navas-Molina JA, Nothias LF, Orchanian SB, Pearson T, Peoples SL, Petras D, Preuss ML, Pruesse E, Rasmussen LB, Rivers A, Robeson, II MS, Rosenthal P, Segata N, Shaffer M, Shiffer A, Sinha R, Song SJ, Spear JR, Swafford AD, Thompson LR, Torres PJ, Trinh P, Tripathi A, Turnbaugh PJ, Ul-Hasan S, van der Hooft JJ, Vargas F, Vázquez-Baeza Y, Vogtmann E, von Hippel M, Walters W, Wan Y, Wang M, Warren J, Weber KC, Williamson CH, Willis AD, Xu ZZ, Zaneveld JR, Zhang Y, Zhu Q, Knight R, Caporaso JG. 2018. QIIME 2: Reproducible, interactive, scalable, and extensible microbiome data science. PeerJ Preprints 6:e27295v2 https://doi.org/10.7287/peerj.preprints.27295v2

译者简介

刘永鑫，博士。2008年毕业于东北农大微生物学专业。2014年中科院遗传发育所获生物信息学博士学位，2016年博士后出站留所工作，任宏基因组学实验室工程师，目前主要研究方向为宏基因组学、数据分析与可重复计算和植物微生物组、QIIME 2项目参与人。发于论文12篇，SCI收录9篇。2017年7月创办“宏基因组”公众号，目前分享宏基因组、扩增子原创文章400+篇，代表博文有《扩增子图表解读、分析流程和统计绘图三部曲》，关注人数3.2万+，累计阅读500万+。

猜你喜欢

系列教程：微生物组入门 Biostar 微生物组  宏基因组

专业技能：学术图表 高分文章 生信宝典 不可或缺的人

一文读懂：宏基因组 寄生虫益处 进化树

必备技能：提问 搜索  Endnote

文献阅读 热心肠 SemanticScholar Geenmedical

扩增子分析：图表解读 分析流程 统计绘图

在线工具：16S预测培养基 生信绘图

科研经验：云笔记  云协作 公众号

编程模板: Shell  R Perl

生物科普:  肠道细菌 人体上的生命 生命大跃进  细胞暗战 人体奥秘  

写在后面

为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内外3000+ 一线科研人员加入。参与讨论，获得专业解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注“姓名-单位-研究方向-职称/年级”。PI请明示身份，另有海内外微生物相关PI群供大佬合作交流。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍末解决群内讨论，问题不私聊，帮助同行。

学习16S扩增子、宏基因组科研思路和分析实战，关注“宏基因组”

点击阅读原文，跳转最新文章目录阅读